/******************************************************************************************/
*  Copyright (c) 2023 HNU/NUIST
*  All rights reserved.
*
*  @Corresponding author :  Qiang Yang(mmmyq@126.com)
*  
*  @author :     De-Gang Chen (degang__chen@126.com)
*  
*  @date	2022-11-2
*
*  @reference : 
*       Junna Zhang, De-Gang Chen, Qiang Yang*, Yiqiao Wang, Dong Liu, Sang-Woon Jeon, and Jun Zhang, "Proximity Ranking-Based Multimodal Differential Evolution", Swarm and Evolutionary Computation, vol. 78, p. 101277, 2023.
*
*  @brief : PRMDE
*  
/******************************************************************************************/

#include "./CEC2013/cec2013.h"
#include "Self_Define_Functions.h"
#include <sys/time.h>
#include <cstdio>
#include <unistd.h>
#include <algorithm>
#include <iostream>
#include <fstream>
#include <stdlib.h>
#include <iomanip>
#include <string>
#include <boost/random.hpp>
#include <boost/random/uniform_int.hpp>
#include <boost/random/cauchy_distribution.hpp>
#include <boost/random/normal_distribution.hpp>
#include <boost/random/uniform_real.hpp>


using namespace std;


// in this program the global best solution is updated after each subgroup finish iterating itself
int main(int argc, char *argv[])
{
    //************Problem parameter************//
    //************Allocate memory for problem parameters************//
    int funToRun[] = {1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20};  //function set
    
    int funNum = 20; //total function num
    

    int population_size_set[] = {100,50,50,50,50,100,1200,350,1200,50,150,200,50,50,200,50,200,300,250,150};


    int i,j,k,fun,run;
    CEC2013 *pFunc = NULL;
    int population_size;
    int dim;
    int Fes;
    int MAX_FES;
    int temp_global_optima_size;
    int global_optima_num;

    int record_count=0;

    int **runs_global_num = new int*[epsilon_set_size];
    int **runs_used_fitness = new int*[epsilon_set_size];

    bool *epsilon_flag = new bool[epsilon_set_size];

    int record_size = 10;

    int *record_fes = new int[record_size];

    int ***runs_global_num_vs_fes = new int**[epsilon_set_size];


    for( i = 0; i < epsilon_set_size; ++i )
    {
        runs_global_num[i] = new int [timesOfRun];
        runs_used_fitness[i] = new int[timesOfRun];

        runs_global_num_vs_fes[i] = new int*[timesOfRun];

        for( j = 0; j < timesOfRun; ++j )
        {
            runs_global_num_vs_fes[i][j] =  new int[record_size];
        }
    }
    //************Problem parameter end************//
    //************Allocate memory for problem parameters end************//

    
    //************DE parameter************//
    //************Allocate memory for DE parameter************//
    int rand_particle_num = 3;      // DE/rand/1 DE/current-to-rand/1
    //int rand_particle_num = 5;    // DE/rand/2 DE/current-to-rand/2
    int *selected_particle_index = new int [rand_particle_num];
    int closest_particle;



    for ( fun =0; fun <funNum; fun ++)
    {
        cout<<"Function "<<funToRun[fun]<<" Begined!"<<endl;

        pFunc = new CEC2013(funToRun[fun]);

        dim = pFunc->get_dimension();

        global_optima_num = pFunc->get_no_goptima();

        MAX_FES = pFunc->get_maxfes();


        double *LBound = new double[dim];
        double *UBound = new double[dim];

        for( i = 0; i < dim; ++i )
        {
            LBound[i] = pFunc->get_lbound(i);
            UBound[i] = pFunc->get_ubound(i);

        }


        for( i = 0; i < record_size; ++i )
        {
            record_fes[i] = (int)((double(i+1))/record_size * MAX_FES);
        }


        population_size = population_size_set[funToRun[fun]-1];
        double **population = new double*[population_size];
        for( i = 0; i < population_size; ++i )
        {
            population[i] = new double[dim];
        }

        double *results = new double [population_size];

        //double *child = new double[dim];
        double **child = new double*[population_size];
        for( i = 0; i < population_size; ++i )
        {
            child[i] = new double[dim];

        }

        double child_fitness;

        Fes = 0;


        // to initialize the population
        boost::mt19937 generator(time(0)*rand());
        boost::uniform_real<> uniform_real_generate_r( 0, 1 );
        boost::variate_generator< boost::mt19937&, boost::uniform_real<> > random_real_num_r( generator, uniform_real_generate_r );

        //************DE parameter end************//
        //************Allocate memory for DE parameter end************//

        //************Algorithm parameter************//
        //************Allocate memory for algorithm parameter************//
        
        double sigma;
        
        //rank_Distance_and_W_and_P
		int **rank_result = new int*[population_size];
	    double **temp_distance = new double*[population_size];
        double **W = new double*[population_size];
        double **P = new double*[population_size];
        for (i = 0; i < population_size; ++i)
		{
			rank_result[i] = new int[population_size];
            temp_distance[i] = new double[population_size];
            W[i] = new double[population_size];
            P[i] = new double[population_size];
		}

        NewType **tempsortdistance = new NewType*[population_size];
        for (i = 0; i < population_size; ++i)
        {
            tempsortdistance[i] = new NewType[population_size];
        }

		//Calculating cumulative fitness
		double **q = new double*[population_size];
		for (i = 0; i < population_size; ++i)
		{
			q[i] = new double[population_size];
		}


        //************DE-algorithm parameter************//
        //************Allocate memory for de-algorithm************//
        //adaptive F-cr
        double *F = new double [population_size];
        double *cr = new double [population_size];
        double *tempresults = new double [population_size];
        double *Fcr = new double [population_size];
        double *Fcr1 = new double [population_size];
        
        NewType *tempsortfitness = new NewType [population_size];
        int *ranking = new int [population_size];

        //************Local search parameter************//
        //************Allocate memory for local search************//
        //PLS
        double *PLSsigma = new double [population_size]; 
        double *PLSP = new double [population_size];
        double *tempresults1 = new double [population_size];
        
        double **PLSxnew = new double*[1];
        for (i = 0; i < 1; ++i)
	    {
		    PLSxnew[i] = new double[dim];
	    }
        double *PLSxnew_fitness = new double[1];

        double *PLSxbetter = new double[dim];
        double PLSxbetter_fitness;
        
        //************Local search parameter end************//
        //************Allocate memory for local search end************//


        //************  timesOfRun  ************//
        for (run = 0; run < timesOfRun; run++)
        {
            cout<<"Running the "<<run<<"th times"<<endl;

            for( i = 0; i < epsilon_set_size; ++i )
            {
                epsilon_flag[i] = false;
            }

            Fes = 0;
            record_count=0;

            // initialize the population
            for( j =0; j < dim; ++j )
            {
                boost::uniform_real<> uniform_real_generate_x( LBound[j], UBound[j] );
                boost::variate_generator< boost::mt19937&, boost::uniform_real<> > random_real_num_x( generator, uniform_real_generate_x );

                for( i =0 ; i < population_size; ++i )
                {
                    population[i][j] = random_real_num_x();
                }
            }

            Popupation_Fitness( population, population_size, Fes, results, pFunc );

		
            while( Fes < MAX_FES )      // Termination condition
            {

                //compute the number of global optima found
                for( i = 0; i < epsilon_set_size; ++i )
                {
                    if( !epsilon_flag[i] )
                    {
                        temp_global_optima_size = Compute_Global_Optima_Found( population,  results, population_size, dim, epsilon_set[i], pFunc );

                        runs_global_num[i][run] = temp_global_optima_size;
                        runs_used_fitness[i][run] = Fes;
                        if( temp_global_optima_size == global_optima_num )
                        {
                            epsilon_flag[i] = true;
                        }
                    }
                }

                if( Fes >= record_fes[record_count] )
                {
                    for( i = 0; i < epsilon_set_size; ++i )
                    {
                        runs_global_num_vs_fes[i][run][record_count] = runs_global_num[i][run];
                    }

                    record_count++;
                }

                // Adaptive calculation of sigma in the weight function
                sigma=(double)(sigmamax-(double)((sigmamax-sigmamin)*(double)(pow((double)Fes/MAX_FES,0.5))));

                // Distance sort and calculate the selection probability
                Ranking_Distance_sort( sigma,population_size, dim, tempsortdistance,rank_result, temp_distance, W,P,q,population);
	
                //adaptiveFcr( population_size, ranking, tempsortfitness, results, cr, F ,tempresults,Fcr);
                //  
                for( i = 0; i < population_size; ++i )
                {
                    tempresults[i] = results[i];
                }
                sort(tempresults,tempresults+population_size,less<double>());
                
                //evlove
                for( i = 0; i < population_size; ++i )
                {
                    // Select individuals to assist the mutation
					Select_Rand_Number( i, population_size, selected_particle_index, rand_particle_num, q,results);

                    // The target vector is different from the base vector, so the individuals serving as guidance vectors for different mutation strategies are different
                    adaptiveFcr1( i,selected_particle_index[0], population_size, ranking, tempsortfitness, results, cr, F ,tempresults,Fcr,Fcr1);

                    // Reproduce the population to the offspring
                    memcpy( child[i], population[i], sizeof(double)*dim );

                    // mutation: DE/rand/1
                    // Note: Theoretically different mutation strategies can be combined here
                    
                    DE_Child(child[i], population[selected_particle_index[0]], population[selected_particle_index[1]], population[selected_particle_index[2]], dim, LBound, UBound, cr[i], F[i] );
                    
                    // mutation: DE/current-to-rand/1
                    //DE_Child1( population[i],child[i], population[selected_particle_index[0]], population[selected_particle_index[1]], population[selected_particle_index[2]], dim, LBound, UBound, cr[i], F[i] );
                    // 

                }

                // CDE select operation
                for( i = 0; i < population_size; ++i )
                {
                    closest_particle = Close_Particle( child[i], population, population_size, dim );

                    child_fitness = Fitness( child[i], Fes, pFunc );

                    if( child_fitness > results[closest_particle] )
                    {
                        results[closest_particle] = child_fitness;
                        memcpy( population[closest_particle], child[i], sizeof(double)*dim );
                    }

                }

                // Adaptive local search
                PLS( population_size, dim,Fes,MAX_FES,ranking,tempsortfitness,tempresults1,results, PLSsigma, PLSxbetter_fitness,PLSP,PLSxbetter,PLSxnew_fitness,PLSxnew,population,LBound, UBound,pFunc );
                           
            }
            //compute the number of global optima found
            for( i = 0; i < epsilon_set_size; ++i )
            {
                if( !epsilon_flag[i] )
                {
                    temp_global_optima_size = Compute_Global_Optima_Found( population, results, population_size, dim, epsilon_set[i], pFunc );

                    runs_global_num[i][run] = temp_global_optima_size;
                    runs_used_fitness[i][run] = MAX_FES;

                }
            }

            for( i = 0; i < epsilon_set_size; ++i )
            {
                runs_global_num_vs_fes[i][run][record_count] = runs_global_num[i][run];
            }


        }

        // Free memory
        for (i = 0; i < population_size; ++i)
		{
			delete[]rank_result[i];
		}
		    delete[]rank_result;
		for (i = 0; i < population_size; ++i)
		{
			delete[]temp_distance[i];
		}
		    delete[]temp_distance;
		for (i = 0; i < population_size; ++i)
		{
			delete[]W[i];
		}
			delete[]W;
		for (i = 0; i < population_size; ++i)
		{
			delete[]P[i];
		}
			delete[]P;

        for (i = 0; i < population_size; ++i)
        {
            delete[]tempsortdistance[i];
        }
        delete []tempsortdistance;
        

		for (i = 0; i < population_size; ++i)
		{
			delete[]q[i];
		}
			delete[]q;
        
        delete[]F;
        delete[]cr;
        delete[]PLSsigma;
        delete[]PLSP;
        delete []tempresults1;
        for (i = 0; i < 1; ++i)
	    {
		    delete[]PLSxnew[i];
	    }
		delete[]PLSxnew;
        delete[]PLSxnew_fitness;
        delete[]PLSxbetter;

        // output
        char fun_name[10];
        char epsilon_counter[10];
        snprintf(fun_name,10,"%d",funToRun[fun]);

        string *Epsilon_Files = new string[epsilon_set_size];
        ofstream *out_Epsilon = new ofstream[epsilon_set_size];

        for( i = 0; i < epsilon_set_size; ++i )
        {
            snprintf(epsilon_counter,10,"%d",i+1);
            Epsilon_Files[i] = "./Results/"+string(epsilon_counter)+"/"+"Final_Optima_Num_And_FES_Epsilon_"+string(epsilon_counter)+"_Fun_"+ string(fun_name)+".txt";
        }

        for( i = 0; i < epsilon_set_size; ++i )
        {
            out_Epsilon[i].open( Epsilon_Files[i].c_str() );
        }

        for( i = 0; i < epsilon_set_size; ++i )
        {
            for( j = 0; j < timesOfRun; ++j )
            {
                out_Epsilon[i]<<runs_global_num[i][j] <<"\t"<<runs_used_fitness[i][j]<<endl;
            }
        }

        for( i = 0; i < epsilon_set_size; ++i )
        {
            out_Epsilon[i].close();
        }

        for( i = 0; i < epsilon_set_size; ++i )
        {
            snprintf(epsilon_counter,10,"%d",i+1);
            Epsilon_Files[i] = "./Results/"+string(epsilon_counter)+"/"+"Optima_Num_VS_FES_Epsilon_"+string(epsilon_counter)+"_Fun_"+ string(fun_name)+".txt";
        }

        for( i = 0; i < epsilon_set_size; ++i )
        {
            out_Epsilon[i].open( Epsilon_Files[i].c_str() );
        }

        for( i = 0; i < epsilon_set_size; ++i )
        {
            for( j = 0; j < timesOfRun; ++j )
            {
                for( k = 0; k < record_size; ++k )
                    out_Epsilon[i]<<runs_global_num_vs_fes[i][j][k] <<"\t";
                out_Epsilon[i]<<endl;
            }
        }

        for( i = 0; i < epsilon_set_size; ++i )
        {
            out_Epsilon[i].close();
        }

        for( i = 0; i < epsilon_set_size; ++i )
        {
            snprintf(epsilon_counter,10,"%d",i+1);
            Epsilon_Files[i] = "./Results/"+string(epsilon_counter)+"/"+"vector_location_"+string(epsilon_counter)+"_Fun_"+ string(fun_name)+".txt";
        }

        for( i = 0; i < epsilon_set_size; ++i )
        {
            out_Epsilon[i].open( Epsilon_Files[i].c_str() );
        }

        for( i = 0; i < epsilon_set_size; ++i )
        {
            for( j = 0; j < population_size; ++j )
            {
                for( k = 0; k < dim; ++k )
		        {
                    out_Epsilon[i]<<population[j][k] <<"\t";
                }
                out_Epsilon[i]<<endl;
            }
        }        

        for( i = 0; i < epsilon_set_size; ++i )
        {
            out_Epsilon[i].close();
        }


        cout<<"Function "<<funToRun[fun]<<" Finished!"<<endl;


        //release the resouces
        for( i =0; i < population_size; ++i )
        {
            delete []population[i];
        }
        delete []population;


        delete []results;
        delete []LBound;
        delete []UBound;
        for( i =0; i < population_size; ++i )
        {
            delete []child[i];
        }
        delete []child;
        delete []Epsilon_Files;
        delete []out_Epsilon;
        
        delete []tempresults;
        delete []Fcr;
        delete []Fcr1;
        
        delete []tempsortfitness;
        delete []ranking;

	}


    for( i = 0; i < epsilon_set_size; ++i )
    {
        delete []runs_global_num[i];
        delete []runs_used_fitness[i];

        for( j = 0; j < timesOfRun; ++j )
        {
            delete []runs_global_num_vs_fes[i][j];
        }
        delete []runs_global_num_vs_fes[i];

    }

    delete []runs_global_num;
    delete []runs_used_fitness;
    delete []epsilon_flag;
    delete []selected_particle_index;

    delete []runs_global_num_vs_fes;

    delete []record_fes;


    delete pFunc;

    return 0;
}